clc;
clear;

%% main主脚本
tic;
startTime = cputime;
pde = process_pde_stokes_t(5);
vector_h = [1/2 1/4 1/8 1/16];
basis_type1 = "P2";
basis_type2 = "P1";


u_Linf = zeros(size(vector_h,2),1);
u_L2 = zeros(size(vector_h,2),1);
u_H1 = zeros(size(vector_h,2),1);
p_Linf = zeros(size(vector_h,2),1);
p_L2 = zeros(size(vector_h,2),1);
p_H1 = zeros(size(vector_h,2),1);

for i = 1:size(vector_h,2)
    h = vector_h(i);
    dt = h^2;
    Nt = ceil((pde.domain(6) - pde.domain(5))/dt);

    mesh = process_mesh(pde, h);
    fs1 = process_fespace(mesh, basis_type1);
    fs2 = process_fespace(mesh, basis_type2);

    % process_showmesh(mesh, fs1);
    % process_showmesh(mesh, fs2);

    M = term_U_V(mesh, fs1);
    M1 = term_P_Q(mesh, fs2);
    A = pde.nu*term_grad_U_grad_V(mesh, fs1);
    U_n = [pde.u1(fs1.node(:,1),fs1.node(:,2),0);
        pde.u2(fs1.node(:,1),fs1.node(:,2),0)];
    s_n = zeros(fs2.Nb,1);
    for m = 0:Nt-1
        t = (m+1)*dt;

        b1 = term_F_V_t(mesh, fs1, pde.f1, pde.f2, t);
        b2 = term_Pold_div_V(mesh, fs1, fs2, s_n);
        Au = 1/dt*M + A;
        bu = 1/dt*M*U_n + b1 + b2;
        for k = 1:size(pde.bdary,2)
            if 0 == pde.bdary(1,k)
                bdn1 = bdary_node(mesh, fs1, k);
                Au(bdn1,:) = 0;
                Au(sub2ind(size(Au), bdn1, bdn1)) = 1;
                bu(bdn1) = pde.u1(fs1.node(bdn1,1),fs1.node(bdn1,2),t);
                bdn2 = bdn1+fs1.Nb;
                Au(bdn2,:) = 0;
                Au(sub2ind(size(Au), bdn2, bdn2)) = 1;
                bu(bdn2) = pde.u2(fs1.node(bdn1,1),fs1.node(bdn1,2),t);
            end
        end
        U_hat = Au\bu;

        Apho = term_grad_P_grad_Q(mesh, fs2);
        bpho = term_div_Uold_Q(mesh, fs1, fs2, U_hat);
        Apho(1,:) = ones(fs2.Nb,1)'*M1;
        bpho(1) = 0;
        pho_n = Apho\bpho;

        As = M1;
        bs = M1*s_n - bpho;
        As(1,:) = ones(fs2.Nb,1)'*M1;
        bs(1) = 0;
        s_n = As\bs;

        b3 = term_grad_Pold_V(mesh, fs1, fs2, pho_n);
        A3 = M;
        b3 = M*U_hat + b3;
        for k = 1:size(pde.bdary,2)
            if 0 == pde.bdary(1,k)
                bdn1 = bdary_node(mesh, fs1, k);
                A3(bdn1,:) = 0;
                A3(sub2ind(size(A3), bdn1, bdn1)) = 1;
                b3(bdn1) = pde.u1(fs1.node(bdn1,1),fs1.node(bdn1,2),t);
                bdn2 = bdn1+fs1.Nb;
                A3(bdn2,:) = 0;
                A3(sub2ind(size(A3), bdn2, bdn2)) = 1;
                b3(bdn2) = pde.u2(fs1.node(bdn1,1),fs1.node(bdn1,2),t);
            end
        end
        U_n = A3\b3;

        P_n = pde.nu*s_n - 1/dt*pho_n;
    end
    U = [U_n; P_n];


    % process_animation(fs1, fs2, U_all, t_all);
    % process_plot(mesh, fs1, fs2, U);
    % process_field(fs1, fs2, U);

    [u_Linf(i,1), p_Linf(i,1)] = process_error_t(mesh, fs1, fs2, pde, U, "Linf");
    [u_L2(i,1), p_L2(i,1)] = process_error_t(mesh, fs1, fs2, pde, U, "L2");
    [u_H1(i,1), p_H1(i,1)] = process_error_t(mesh, fs1, fs2, pde, U, "H1");

    fprintf("已完成\t%d/%d\n", i, size(vector_h,2));
end

process_rate(u_Linf, vector_h);
process_rate(u_L2, vector_h);
process_rate(u_H1, vector_h);
process_rate(p_Linf, vector_h);
process_rate(p_L2, vector_h);
process_rate(p_H1, vector_h);

TOCtime = toc;
CPUtime = cputime - startTime;
disp(["TOC time(s)", num2str(TOCtime)]);
disp(["CPU time(s)", num2str(CPUtime)]);

